Signal Decay¶
Important
This chapter should differentiate itself from TE_Dependence by focusing on the application of (1) to raw signal, rather than decompositions.
Important
The general flow of this chapter should be:
Explain the monoexponential decay equation.
Mention other models of signal decay, including biexponential models, and why they generally aren’t used on multi-echo fMRI data.
Show signal decay in real data.
Start with single-volume brain images across increasing echo times.
Next, show a scatter plot of the values for a single voxel over time, separated by echo time, to give a better sense of the scale of the decay.
Line plots of the different echo times for a single voxel to show that decay is not a simple scalar.
Show a simulated, static decay curve.
Show the difference in decay curves for “active” and “inactive” states in the same voxel (i.e., different T2* values).
Show how single-echo data fluctuate over time, given an underlying signal.
Add a changing decay curve underneath the single-echo data point.
Overlay points for the changing T2* and S0 values as well.
Plot the decay curves of the same relative fluctuations, but dependent solely on S0 or T2* changes.
Overlay single-echo data points on the fluctuating decay curves, to show how the two types of fluctuations cannot be dissociated using single-echo data.
Plot scaled T2* against the single-echo data in a line plot. The time series will be very similar, but slightly different.
Plot scaled S0 against the single-echo data. The time series will effectively be the same.
Point out that we don’t know the underlying T2* or S0 fluctuations, so this just shows how very similar signal can come from T2* or S0.
Overlay multi-echo data points on the fluctuating decay curves, to show how the differences in the patterns are much clearer.
Discuss the effect of noise on the multi-echo signal. Perhaps plot a single voxel’s real values?
Point to the TE Dependence chapter.
In this chapter, we illustrate how BOLD signal decays and how we can glean useful information from that decay rate.
To do this, we will use a combination of simulated and real multi-echo data.
import os
import matplotlib.pyplot as plt
import nibabel as nib
import numpy as np
import pandas as pd
import seaborn as sns
from IPython import display
from matplotlib.animation import FuncAnimation
from myst_nb import glue
from nilearn import image, masking, plotting
from nilearn.glm import first_level
from repo2data.repo2data import Repo2Data
from scipy import signal
from book_utils import predict_bold_signal
sns.set_style("whitegrid")
plt.rcParams.update(
{"text.usetex": True, "font.family": "sans-serif", "font.sans-serif": ["Helvetica"]}
)
# Install the data if running locally, or point to cached data if running on neurolibre
DATA_REQ_FILE = os.path.join("../binder/data_requirement.json")
# Download data
repo2data = Repo2Data(DATA_REQ_FILE)
data_path = repo2data.install()
data_path = os.path.abspath(os.path.join(data_path[0], "data"))
out_dir = os.path.join(data_path, "signal-decay")
os.makedirs(out_dir, exist_ok=True)
---- repo2data starting ----
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/repo2data
Config from file :
../binder/data_requirement.json
Destination:
./../data/multi-echo-data-analysis
Info : ./../data/multi-echo-data-analysis already downloaded
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/nilearn/datasets/__init__.py:96: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
"Numpy arrays.", FutureWarning)
func_dir = os.path.join(data_path, "sub-04570/func/")
data_files = [
os.path.join(
func_dir,
"sub-04570_task-rest_echo-1_space-scanner_desc-partialPreproc_bold.nii.gz",
),
os.path.join(
func_dir,
"sub-04570_task-rest_echo-2_space-scanner_desc-partialPreproc_bold.nii.gz",
),
os.path.join(
func_dir,
"sub-04570_task-rest_echo-3_space-scanner_desc-partialPreproc_bold.nii.gz",
),
os.path.join(
func_dir,
"sub-04570_task-rest_echo-4_space-scanner_desc-partialPreproc_bold.nii.gz",
),
]
ted_dir = os.path.join(data_path, "tedana")
ECHO_TIMES = [12, 28, 44, 60]
n_echoes = len(ECHO_TIMES)
pal = sns.color_palette("cubehelix", n_echoes)
img = image.index_img(data_files[0], 0)
data = img.get_fdata()
vmax = np.max(data)
idx = np.vstack(np.where(data > 1000))
min_x, min_y, min_z = np.min(idx, axis=1)
max_x, max_y, max_z = np.max(idx, axis=1)
imgs = []
for f in data_files:
img = image.index_img(f, 0)
data = img.get_fdata()
data = data[min_x:max_x, min_y:max_y, min_z:max_z]
img = nib.Nifti1Image(data, img.affine, img.header)
imgs.append(img)
# Create binary mask from adaptive mask
adaptive_mask_file = os.path.join(
ted_dir, "sub-04570_task-rest_space-scanner_desc-adaptiveGoodSignal_mask.nii.gz"
)
mask = image.math_img("img >= 3", img=adaptive_mask_file)
# Component parameter estimates
betas_file = os.path.join(
ted_dir, "sub-04570_task-rest_space-scanner_desc-ICA_components.nii.gz"
)
beta_maps = masking.apply_mask(betas_file, mask)
# Component table
comp_tbl = pd.read_table(
os.path.join(ted_dir, "sub-04570_task-rest_space-scanner_desc-tedana_metrics.tsv"),
index_col="Component",
)
# Get voxel index for voxel most related to component with highest kappa value
acc_comp_tbl = comp_tbl.loc[comp_tbl["classification"] == "accepted"]
high_kappa_comp = acc_comp_tbl.sort_values(by="kappa", ascending=False).index.values[0]
high_kappa_comp_val = int(high_kappa_comp.split("_")[1])
voxel_idx = np.where(
beta_maps[high_kappa_comp_val, :] == np.max(beta_maps[high_kappa_comp_val, :])
)[0][0]
# Extract the time series for the high-kappa voxel
data = [masking.apply_mask(f, mask) for f in data_files]
ts = [d[:, voxel_idx] for d in data]
ts_1d = np.hstack(ts)
n_trs = ts[0].shape[0]
Signal decays as echo time increases¶
To start, let us look at how BOLD signal changes with increasing echo time in real data. The data we use is a single resting-state run, with echo times of 12, 28, 44, and 60 milliseconds.
plt.style.use("dark_background")
fig, axes = plt.subplots(figsize=(26, 4), ncols=len(data_files))
for i_echo, img in enumerate(imgs):
te = ECHO_TIMES[i_echo]
if i_echo == 0:
data = img.get_fdata()
vmax = np.max(data)
plotting.plot_epi(
img,
cut_coords=[0],
display_mode="z",
annotate=False,
draw_cross=False,
axes=axes[i_echo],
figure=fig,
vmax=vmax,
cmap="gray",
)
axes[i_echo].set_title("TE={}ms".format(te), fontsize=20, pad=0)
glue("fig_brain_decay", fig, display=False)
# Reset the style
plt.style.use("default")
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
253 command, cwd=cwd if cwd is not None else self.texcache,
--> 254 stderr=subprocess.STDOUT)
255 except FileNotFoundError as exc:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in check_output(timeout, *popenargs, **kwargs)
410 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411 **kwargs).stdout
412
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
487
--> 488 with Popen(*popenargs, **kwargs) as process:
489 try:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)
799 errread, errwrite,
--> 800 restore_signals, start_new_session)
801 except:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
1550 err_msg += ': ' + repr(err_filename)
-> 1551 raise child_exception_type(errno_num, err_msg, err_filename)
1552 raise child_exception_type(err_msg)
FileNotFoundError: [Errno 2] No such file or directory: 'latex': 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
339 pass
340 else:
--> 341 return printer(obj)
342 # Finally look for special method names
343 method = get_real_method(obj, self.print_method)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
149 FigureCanvasBase(fig)
150
--> 151 fig.canvas.print_figure(bytes_io, **kw)
152 data = bytes_io.getvalue()
153 if fmt == 'svg':
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2228 else suppress())
2229 with ctx:
-> 2230 self.figure.draw(renderer)
2231
2232 if bbox_inches:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
72 @wraps(draw)
73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74 result = draw(artist, renderer, *args, **kwargs)
75 if renderer._rasterizing:
76 renderer.stop_rasterizing()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
2779 self.patch.draw(renderer)
2780 mimage._draw_list_compositing_images(
-> 2781 renderer, self, artists, self.suppressComposite)
2782
2783 for sfig in self.subfigs:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/_api/deprecation.py in wrapper(*inner_args, **inner_kwargs)
429 else deprecation_addendum,
430 **kwargs)
--> 431 return func(*inner_args, **inner_kwargs)
432
433 return wrapper
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
2879 artists.remove(spine)
2880
-> 2881 self._update_title_position(renderer)
2882
2883 if not self.axison or inframe:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in _update_title_position(self, renderer)
2810 if (ax.xaxis.get_ticks_position() in ['top', 'unknown']
2811 or ax.xaxis.get_label_position() == 'top'):
-> 2812 bb = ax.xaxis.get_tightbbox(renderer)
2813 else:
2814 bb = ax.get_window_extent(renderer)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in get_tightbbox(self, renderer, for_layout_only)
1081 ticks_to_draw = self._update_ticks()
1082
-> 1083 self._update_label_position(renderer)
1084
1085 # go back to just this axis's tick labels
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _update_label_position(self, renderer)
2078 # get bounding boxes for this axis and any siblings
2079 # that have been set by `fig.align_xlabels()`
-> 2080 bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
2081
2082 x, y = self.label.get_position()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_boxes_siblings(self, renderer)
1866 axis = ax._get_axis_map()[axis_name]
1867 ticks_to_draw = axis._update_ticks()
-> 1868 tlb, tlb2 = axis._get_tick_bboxes(ticks_to_draw, renderer)
1869 bboxes.extend(tlb)
1870 bboxes2.extend(tlb2)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_bboxes(self, ticks, renderer)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in <listcomp>(.0)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in get_window_extent(self, renderer, dpi)
901
902 with cbook._setattr_cm(self.figure, dpi=dpi):
--> 903 bbox, info, descent = self._get_layout(self._renderer)
904 x, y = self.get_unitless_position()
905 x, y = self.get_transform().transform((x, y))
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in _get_layout(self, renderer)
306 _, lp_h, lp_d = renderer.get_text_width_height_descent(
307 "lp", self._fontproperties,
--> 308 ismath="TeX" if self.get_usetex() else False)
309 min_dy = (lp_h - lp_d) * self._linespacing
310
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
228 fontsize = prop.get_size_in_points()
229 w, h, d = texmanager.get_text_width_height_descent(
--> 230 s, fontsize, renderer=self)
231 return w, h, d
232
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
397 else:
398 # use dviread.
--> 399 dvifile = self.make_dvi(tex, fontsize)
400 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
401 page, = dvi
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in make_dvi(self, tex, fontsize)
291 self._run_checked_subprocess(
292 ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 293 texfile], tex, cwd=tmpdir)
294 (Path(tmpdir) / Path(dvifile).name).replace(dvifile)
295 return dvifile
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
256 raise RuntimeError(
257 'Failed to process string with tex because {} could not be '
--> 258 'found'.format(command[0])) from exc
259 except subprocess.CalledProcessError as exc:
260 raise RuntimeError(
RuntimeError: Failed to process string with tex because latex could not be found
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
253 command, cwd=cwd if cwd is not None else self.texcache,
--> 254 stderr=subprocess.STDOUT)
255 except FileNotFoundError as exc:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in check_output(timeout, *popenargs, **kwargs)
410 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411 **kwargs).stdout
412
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
487
--> 488 with Popen(*popenargs, **kwargs) as process:
489 try:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)
799 errread, errwrite,
--> 800 restore_signals, start_new_session)
801 except:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
1550 err_msg += ': ' + repr(err_filename)
-> 1551 raise child_exception_type(errno_num, err_msg, err_filename)
1552 raise child_exception_type(err_msg)
FileNotFoundError: [Errno 2] No such file or directory: 'latex': 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
339 pass
340 else:
--> 341 return printer(obj)
342 # Finally look for special method names
343 method = get_real_method(obj, self.print_method)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
149 FigureCanvasBase(fig)
150
--> 151 fig.canvas.print_figure(bytes_io, **kw)
152 data = bytes_io.getvalue()
153 if fmt == 'svg':
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2228 else suppress())
2229 with ctx:
-> 2230 self.figure.draw(renderer)
2231
2232 if bbox_inches:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
72 @wraps(draw)
73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74 result = draw(artist, renderer, *args, **kwargs)
75 if renderer._rasterizing:
76 renderer.stop_rasterizing()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
2779 self.patch.draw(renderer)
2780 mimage._draw_list_compositing_images(
-> 2781 renderer, self, artists, self.suppressComposite)
2782
2783 for sfig in self.subfigs:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/_api/deprecation.py in wrapper(*inner_args, **inner_kwargs)
429 else deprecation_addendum,
430 **kwargs)
--> 431 return func(*inner_args, **inner_kwargs)
432
433 return wrapper
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
2879 artists.remove(spine)
2880
-> 2881 self._update_title_position(renderer)
2882
2883 if not self.axison or inframe:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in _update_title_position(self, renderer)
2810 if (ax.xaxis.get_ticks_position() in ['top', 'unknown']
2811 or ax.xaxis.get_label_position() == 'top'):
-> 2812 bb = ax.xaxis.get_tightbbox(renderer)
2813 else:
2814 bb = ax.get_window_extent(renderer)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in get_tightbbox(self, renderer, for_layout_only)
1081 ticks_to_draw = self._update_ticks()
1082
-> 1083 self._update_label_position(renderer)
1084
1085 # go back to just this axis's tick labels
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _update_label_position(self, renderer)
2078 # get bounding boxes for this axis and any siblings
2079 # that have been set by `fig.align_xlabels()`
-> 2080 bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
2081
2082 x, y = self.label.get_position()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_boxes_siblings(self, renderer)
1866 axis = ax._get_axis_map()[axis_name]
1867 ticks_to_draw = axis._update_ticks()
-> 1868 tlb, tlb2 = axis._get_tick_bboxes(ticks_to_draw, renderer)
1869 bboxes.extend(tlb)
1870 bboxes2.extend(tlb2)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_bboxes(self, ticks, renderer)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in <listcomp>(.0)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in get_window_extent(self, renderer, dpi)
901
902 with cbook._setattr_cm(self.figure, dpi=dpi):
--> 903 bbox, info, descent = self._get_layout(self._renderer)
904 x, y = self.get_unitless_position()
905 x, y = self.get_transform().transform((x, y))
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in _get_layout(self, renderer)
306 _, lp_h, lp_d = renderer.get_text_width_height_descent(
307 "lp", self._fontproperties,
--> 308 ismath="TeX" if self.get_usetex() else False)
309 min_dy = (lp_h - lp_d) * self._linespacing
310
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
228 fontsize = prop.get_size_in_points()
229 w, h, d = texmanager.get_text_width_height_descent(
--> 230 s, fontsize, renderer=self)
231 return w, h, d
232
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
397 else:
398 # use dviread.
--> 399 dvifile = self.make_dvi(tex, fontsize)
400 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
401 page, = dvi
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in make_dvi(self, tex, fontsize)
291 self._run_checked_subprocess(
292 ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 293 texfile], tex, cwd=tmpdir)
294 (Path(tmpdir) / Path(dvifile).name).replace(dvifile)
295 return dvifile
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
256 raise RuntimeError(
257 'Failed to process string with tex because {} could not be '
--> 258 'found'.format(command[0])) from exc
259 except subprocess.CalledProcessError as exc:
260 raise RuntimeError(
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 1872x288 with 8 Axes>
<Figure size 1872x288 with 8 Axes>
Fig. 1 Signal decay in the brain.¶
Echo-specific data and echo time¶
fig, ax = plt.subplots(figsize=(14, 6))
values = [i[0] for i in ts]
for i_echo in range(n_echoes):
rep_echo_times = np.ones(n_trs) * ECHO_TIMES[i_echo]
ax.scatter(rep_echo_times, ts[i_echo], alpha=0.05, color=pal[i_echo])
ax.set_ylabel("BOLD signal", fontsize=16)
ax.set_xlabel("Echo Time (ms)", fontsize=16)
ax.set_xticks(ECHO_TIMES)
ax.tick_params(axis="both", which="major", labelsize=14)
ax.set_xlim(0, 70)
ax.set_ylim(0, 3000)
fig.tight_layout()
glue("fig_echo_scatter", fig, display=False)
Fig. 2 Scatter plot of voxel’s values by echo time.¶
The signal itself changes with echo time as well¶
While the overall scale of the signal decreases with echo time, the actual pattern of the signal changes as well.
fig, axes = plt.subplots(n_echoes, sharex=True, sharey=False, figsize=(14, 6))
for i_echo in range(n_echoes):
axes[i_echo].plot(ts[i_echo], color=pal[i_echo])
axes[i_echo].set_ylabel(
f"{ECHO_TIMES[i_echo]}ms", rotation=0, va="center", ha="right", fontsize=14
)
axes[i_echo].set_yticks([])
axes[i_echo].set_xticks([])
axes[-1].set_xlabel("Time", fontsize=16)
axes[-1].set_xlim(0, len(ts[i_echo]) - 1)
fig.tight_layout()
glue("fig_echo_timeseries", fig, display=False)
Fig. 3 Time series from a voxel for each echo.¶
The impact of \(T_{2}^{*}\) and \(S_{0}\) fluctuations on BOLD signal¶
In this section, we investigate the two factors driving changes in the signal decay pattern: \(T_{2}^{*}\) and \(S_{0}\).
# Simulate data
SINGLEECHO_TE = np.array([30])
# For a nice, smooth curve
FULLCURVE_TES = np.arange(0, 101, 1)
n_echoes = len(FULLCURVE_TES)
pal = sns.color_palette("cubehelix", 8)
Simulate \(T_{2}^{*}\) and \(S_{0}\) fluctuations¶
# Simulate data
# We'll convolve with HRF just for smoothness
hrf = first_level.spm_hrf(1, oversampling=1)
N_VOLS = 21
SCALING_FRACTION = 0.1 # used to scale standard deviation
MEAN_T2S = 35
MEAN_S0 = 16000
# simulate the T2*/S0 time series
# The original time series will be a random time series from a normal distribution,
# convolved with the HRF
ts = np.random.normal(loc=0, scale=1, size=(N_VOLS + 20,))
ts = signal.convolve(ts, hrf)[20 : N_VOLS + 20]
ts *= SCALING_FRACTION / np.std(ts)
ts -= ts[0]
t2s_ts = (ts * MEAN_T2S) + MEAN_T2S
s0_ts = (ts * MEAN_S0) + MEAN_S0
Plot BOLD signal decay for a standard single-echo scan¶
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [30])
fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
)
ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
glue("fig_signal_decay_single-echo", fig, display=False)
Fig. 4 BOLD signal decay for a standard single-echo scan¶
Plot BOLD signal decay and BOLD contrast¶
fullcurve_signal_active = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [40])
fullcurve_signal_inactive = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [20])
fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
FULLCURVE_TES,
fullcurve_signal_active[:, 0],
alpha=0.5,
color="red",
label="High Activity",
)
ax.plot(
FULLCURVE_TES,
fullcurve_signal_inactive[:, 0],
alpha=0.5,
color="blue",
label="Low Activity",
)
ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.legend(fontsize=20)
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
glue("fig_signal_decay_contrast", fig, display=False)
Fig. 5 BOLD signal decay and BOLD contrast¶
Plot single-echo data resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how fMRI data fluctuates over time.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
# gif code based on https://www.geeksforgeeks.org/create-an-animated-gif-using-python-matplotlib/
fig, axes = plt.subplots(
nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]}
)
# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(ECHO_TIMES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
ax0_line_plot = axes[0].plot(singleecho_signal[0, :], color="black", zorder=0)
ax0_scatter_plot = axes[0].scatter(
0,
singleecho_signal[:, 0],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, 0],
color="orange",
s=150,
alpha=1.0,
label="Single-Echo Signal",
zorder=1,
)
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack(
(
frame,
singleecho_signal[:, frame],
)
)
)
ax1_scatter_plot.set_offsets(
np.column_stack(
(
SINGLEECHO_TE,
singleecho_signal[:, frame],
)
)
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_single-echo", html, display=False)
Fig. 6 Single-echo data resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how single-echo data is a sample from a signal decay curve.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
# gif code based on https://www.geeksforgeeks.org/create-an-animated-gif-using-python-matplotlib/
fig, axes = plt.subplots(
nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]}
)
# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(ECHO_TIMES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
ax0_line_plot = axes[0].plot(singleecho_signal[0, :], color="black", zorder=0)
ax0_scatter_plot = axes[0].scatter(
0,
singleecho_signal[:, 0],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_line_plot = axes[1].plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
zorder=0,
)[0]
ax1_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, 0],
color="orange",
s=150,
alpha=1.0,
label="Single-Echo Signal",
zorder=1,
)
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack(
(
frame,
singleecho_signal[:, frame],
)
)
)
ax1_line_plot.set_data(
(
FULLCURVE_TES,
fullcurve_signal[:, frame],
)
)
ax1_scatter_plot.set_offsets(
np.column_stack(
(
SINGLEECHO_TE,
singleecho_signal[:, frame],
)
)
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay", html, display=False)
Fig. 7 Single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how changes in fMRI data can be driven by both S0 and T2* fluctuations.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
fig, axes = plt.subplots(
nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]}
)
t2s_value = predict_bold_signal(
np.array([t2s_ts[0]]), np.array([s0_ts[0]]), np.array([t2s_ts[0]])
)[0]
ax0_line_plot = axes[0].plot(singleecho_signal[0, :], color="black", zorder=0)
ax0_scatter_plot = axes[0].scatter(
0,
singleecho_signal[:, 0],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_line_plot = axes[1].plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
zorder=0,
)[0]
ax1_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, 0],
color="orange",
s=150,
alpha=1.0,
label="Single-Echo Signal",
zorder=1,
)
ax1_t2s_scatter_plot = axes[1].scatter(
t2s_ts[0],
t2s_value,
marker="*",
color="blue",
s=500,
alpha=0.5,
label="$T_{2}^{*}$",
zorder=1,
)
ax1_s0_scatter_plot = axes[1].scatter(
0,
s0_ts[0],
marker="*",
color="red",
s=500,
alpha=0.5,
label="$S_{0}$",
clip_on=False,
zorder=2,
)
# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(ECHO_TIMES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
t2s_value = predict_bold_signal(
np.array([t2s_ts[frame]]), np.array([s0_ts[frame]]), np.array([t2s_ts[frame]])
)[0]
ax0_scatter_plot.set_offsets(
np.column_stack(
(
frame,
singleecho_signal[:, frame],
)
)
)
ax1_line_plot.set_data(
(
FULLCURVE_TES,
fullcurve_signal[:, frame],
)
)
ax1_scatter_plot.set_offsets(
np.column_stack(
(
SINGLEECHO_TE,
singleecho_signal[:, frame],
)
)
)
ax1_t2s_scatter_plot.set_offsets(
np.column_stack(
(
t2s_ts[frame],
t2s_value,
)
)
)
ax1_s0_scatter_plot.set_offsets(
np.column_stack(
(
0,
s0_ts[frame],
)
)
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay2", html, display=False)
Fig. 8 Single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how fluctuations in S0 and T2* produce different patterns in the full signal decay curves.
s0based_fullcurve_signal = predict_bold_signal(
FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S)
)
t2sbased_fullcurve_signal = predict_bold_signal(
FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts
)
unperturbed_fullcurve_signal = np.mean(
(s0based_fullcurve_signal[:, 0], t2sbased_fullcurve_signal[:, 0]),
axis=0,
)
fig, axes = plt.subplots(
nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]}
)
t2s_value = predict_bold_signal(
np.array([t2s_ts[0]]), np.array([s0_ts[0]]), np.array([t2s_ts[0]])
)[0]
ax0_line_plot = axes[0].plot(ts, color="black", zorder=0)
ax0_scatter_plot = axes[0].scatter(
0,
ts[0],
color="purple",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_orig_line_plot = axes[1].plot(
FULLCURVE_TES,
unperturbed_fullcurve_signal,
alpha=0.5,
color="black",
label="Original Average Decay Curve",
zorder=-1,
)[0]
ax1_s0_line_plot = axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, 0],
alpha=0.5,
color="red",
label="$S_0$-Driven Decay Curve",
)[0]
ax1_t2s_line_plot = axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, 0],
alpha=0.5,
color="blue",
label="$T_{2}^{*}$-Driven Decay Curve",
)[0]
# Set constant params for figure
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(ECHO_TIMES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack(
(
frame,
ts[frame],
)
)
)
ax1_t2s_line_plot.set_data(
(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, frame],
)
)
ax1_s0_line_plot.set_data(
(
FULLCURVE_TES,
s0based_fullcurve_signal[:, frame],
)
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay3", html, display=False)
Fig. 9 \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data¶
This shows how single-echo data, on its own, cannot distinguish between S0 and T2* fluctuations.
s0based_fullcurve_signal = predict_bold_signal(
FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S)
)
t2sbased_fullcurve_signal = predict_bold_signal(
FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts
)
s0based_singleecho_signal = s0based_fullcurve_signal[SINGLEECHO_TE, :]
t2sbased_singleecho_signal = t2sbased_fullcurve_signal[SINGLEECHO_TE, :]
fig, axes = plt.subplots(
nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]}
)
ax0_line_plot = axes[0].plot(ts, color="black", zorder=0)
ax0_scatter_plot = axes[0].scatter(
0,
ts[0],
color="purple",
s=150,
zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, 0],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
zorder=0,
)[0]
ax1_t2s_line_plot = axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, 0],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
zorder=1,
)[0]
ax1_s0_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
s0based_singleecho_signal[:, 0],
color="red",
s=150,
alpha=1.0,
label="$S_0$-Driven Single-Echo Signal",
zorder=2,
)
ax1_t2s_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
t2sbased_singleecho_signal[:, 0],
color="blue",
s=150,
alpha=1.0,
label="$T_{2}^{*}$-Driven Single-Echo Signal",
zorder=3,
)
# Set constant params for figure
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(ECHO_TIMES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack(
(
frame,
ts[frame],
)
)
)
ax1_s0_line_plot.set_data(
(
FULLCURVE_TES,
s0based_fullcurve_signal[:, frame],
)
)
ax1_t2s_line_plot.set_data(
(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, frame],
)
)
ax1_s0_scatter_plot.set_offsets(
np.column_stack(
(
SINGLEECHO_TE,
s0based_singleecho_signal[:, frame],
)
)
)
ax1_t2s_scatter_plot.set_offsets(
np.column_stack(
(
SINGLEECHO_TE,
t2sbased_singleecho_signal[:, frame],
)
)
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay4", html, display=False)
Fig. 10 \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data¶
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data¶
This shows how S0 and T2* fluctuations produce different patterns in multi-echo data.
s0based_fullcurve_signal = predict_bold_signal(
FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S)
)
t2sbased_fullcurve_signal = predict_bold_signal(
FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts
)
s0based_multiecho_signal = s0based_fullcurve_signal[ECHO_TIMES, :]
t2sbased_multiecho_signal = t2sbased_fullcurve_signal[ECHO_TIMES, :]
fig, axes = plt.subplots(
nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]}
)
ax0_line_plot = axes[0].plot(
ts,
color="black",
zorder=0,
)[0]
ax0_scatter_plot = axes[0].scatter(
0,
ts[0],
color="purple",
s=150,
zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, 0],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
zorder=0,
)[0]
ax1_t2s_line_plot = axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, 0],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
zorder=1,
)[0]
ax1_s0_scatter_plot = axes[1].scatter(
ECHO_TIMES,
s0based_multiecho_signal[:, 0],
color="red",
s=150,
alpha=1.0,
label="$S_0$-Driven Multi-Echo Signal",
zorder=2,
)
ax1_t2s_scatter_plot = axes[1].scatter(
ECHO_TIMES,
t2sbased_multiecho_signal[:, 0],
color="blue",
s=150,
alpha=1.0,
label="$T_{2}^{*}$-Driven Multi-Echo Signal",
zorder=3,
)
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(ECHO_TIMES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack(
(
frame,
ts[frame],
)
)
)
ax1_s0_line_plot.set_data(
(
FULLCURVE_TES,
s0based_fullcurve_signal[:, frame],
)
)
ax1_t2s_line_plot.set_data(
(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, frame],
)
)
ax1_s0_scatter_plot.set_offsets(
np.column_stack(
(
ECHO_TIMES,
s0based_multiecho_signal[:, frame],
)
)
)
ax1_t2s_scatter_plot.set_offsets(
np.column_stack(
(
ECHO_TIMES,
t2sbased_multiecho_signal[:, frame],
)
)
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay5", html, display=False)
Fig. 11 \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data¶
Plot \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)¶
When the BOLD signal is driven entirely by T2* fluctuations, the BOLD signal is very similar to a scaled version of those T2* fluctuations, but there are small differences.
fig, ax = plt.subplots(figsize=(16, 4))
scalar = np.linalg.lstsq(
t2s_ts[:, None], t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :], rcond=None
)[0]
ax.plot(
t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
color="orange",
label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
t2s_ts * scalar,
label="$T_{2}^{*}$ (scaled to signal)",
linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
glue("fig_t2s_bold_single-echo", fig, display=False)
Fig. 12 \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)¶
Plot \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)¶
When the BOLD signal is driven entirely by S0 fluctuations, the BOLD signal ends up being a scaled version of the S0 fluctuations.
fig, ax = plt.subplots(figsize=(16, 4))
scalar = np.linalg.lstsq(
s0_ts[:, None], s0based_fullcurve_signal[SINGLEECHO_TE[0], :], rcond=None
)[0]
ax.plot(
s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
color="orange",
label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
s0_ts * scalar,
label="S0 (scaled to signal)",
linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
glue("fig_s0_bold_single-echo", fig, display=False)
Fig. 13 \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)¶